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Abstract 

We analyze a numerical instability that occurs in the well-known split-step 
Fourier method on the background of a soliton. This instability is found to be 
very sensitive to small changes of the parameters of both the numerical grid and 
the soliton, unlike the instability of most finite-difference schemes. Moreover, 
the principle of "frozen coefficients" , in which variable coefficients are treated as 
"locally constant" for the purpose of stability analysis, is strongly violated for 
the instability of the split-step method on the soliton background. Our analysis 
explains all these features. It is enabled by the fact that the period of oscillations 
of the unstable Fourier modes is much smaller than the width of the soliton. 
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1 Introduction 



The split-step Fourier method is widely used in numerical simulations of nonlinear 
wave equations. In particular, it is the mainstream method in nonlinear optics, where 
the fundamental equation describing propagation of an electromagnetic pulse or beam 
is the nonlinear Schrodinger equation 

iuz - /3utt + ^ulul"^ = , u{t,0) = uo{t). (1.1) 

In this paper, we use notations adopted in fiber optics [T], whereby u is proportional to 
the complex envelope of the electric field, z is the propagation distance along the fiber, 
and t is the time in the reference frame moving with the pulse. (The subscripts denote 
partial differentiation.) Thus, z and t are the evolution and spatial variables, respec- 
tively. Also, in Eq. (11. ip . /3 and 7 are proportional to the group velocity dispersion and 
nonlinear refractive index of the optical fiber, respectively. Although these real-valued 
constants can be scaled out of the equation by an appropriate nondimensionalization, 
we will keep them in our analysis to distinguish contributions of the dispersive and 
nonlinear terms. Solitons with u{z, \t\ — 00) — exist in (II. ip for (3'-f < 0. 

The idea of the split-step method is the following. Equation (11. ip can be solved 
analytically when only one of the dispersive and nonlinear terms is nonzero. Then, the 
approximate numerical solution of (II. ip can be obtained in a sequence of steps which 
alternatingly account either for the dispersion or for the nonlinearity: 



for n from 1 to rimax do: 

Mnonlm(t) = Unit) CXp (i7|M„(t) ^ Az) 
Un+1 = exp(2/3a;^ A2;) J'[unonlin(t)] 

end 



1.2) 



Here A2; is the step size along the evolution variable, and the subscript n now denotes 
the value at the nth step0. In algorithm (II. 2p . the evolution due to the dispersive term 
is computed using the Fourier transform T and its inverse where, e.g.: 



7\u\{u)= \ u{t)e-'^'dt. (1.3) 
J —00 

For equations of the form (II. ip . such an algorithm was originally applied in [2] - |1] 
and later comprehensively studied in [5] . Implementations where the dispersive term is 
computed using finite-difference discretization in t can be used as well, but the Fourier- 
based implementation is preferred in optics, not only because of its higher accuracy 
(exponential versus algebraic in At for smooth pulses), but also because it allows one 
to easily handle more complicated dispersive terms than in (II. ip . In this paper, we 
specifically focus on the Fourier-based implementation (II. 2p of the split-step method. 



^Using subscripts to denote iteration steps, as in (|1.2p . and partial differentiation, as in (jl.ip . will 
not lead to confusion. 
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Since the split-step method is exphcit, it can only be conditionally stable. However, 
the standard von Neumann stability analysis on the background of the trivial solution 
u{t) = does not reveal any instability. Weideman and Herbst [6] were the first to 
rigorously show that the split-step method is indeed conditionally stable, with the in- 
stability being able to develop over the background of a finite- amplitude monochromatic 
wave 

Ucw = A exp{iKz -i^c^t), = + 7|A|^ A = const. (1.4) 

A few years later, Matera et al [7] obtained a similar result by a more heuristic method. 
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Figure 1: Schematics of the numerical instability increment of method (11. 2p over the 
background of a monocromatic wave (ll.4p as a function of frequency u. Only the part 
for > is shown, since the graph is symmetric about u = 0. The width of the 
instability peak near Ukn, where k = 1,2, . . ., equals I'jA'^ / {f3ukn)\ (see Appendix 3). 
Note that the case corresponding to /3 < is shown; in the case of /3 > 0, the peaks 
will occur on the opposite sides of cOt, and U2tt- 



The numerical instability increment of the split-step method of Eq. (II. ip with the 
background solution u = Ucw (11. 4p . found by Weideman and Herbst, is schematically 
shown in Fig. [1] for a particular case of Qcw = 0. (The formulae describing the location 
and growth rate of unstable modes are given in Appendix 3.) Here the notations are 
the following. The numerical solution of (II. ip is assumed to have the form 



Ur. 



Un = U + Un, 



Aexp{\zn — iut), 



\Un\ <^ \u\ 



Zr,. 



nAz, u 



(1.5a) 
^1.5b) 



where —T/2 <t<T/2 and the limits for the integer index I are determined by the 
number of grid points. In Fig. [H 



TT 



i.e. \I3\ujIAz = 71, 



1.6) 
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and, similarly, |/3|co'27r'^-^ — ^tt. Instability peaks of the same height exist also near 
iwsTT, as long as cjmax = ^r/At, where At is the mesh size along t, extends beyond cjs^; 
etc. 

Note that the instability depicted in Fig. [1] is mild, in the sense that the numerical 
error grows by a factor exp (0(1)) over propagation distances of order 0(1). Further- 
more, this instability remains mild even when Az significantly exceeds its threshold 

TT 

A2;thresh ~ p (1-7) 

(see (11. 6p ). This should be contrasted with the behavior of finite-difference explicit 
schemes, e.g., the Runge-Kutta scheme for the Heat equation, where the instability be- 
comes strong (i.e. the numerical error grows by a factor of exp (0(1)) over propagation 
distances of only 0{Az)) whenever the step size along the evolution variable exceeds 
the instability threshold. Another remarkable property of the instability in Fig. [1] is 
that it occurs only in a finite (and rather narrow, of width 0{'yA'^ /(\(3\u.j^) = 0{V Az)) 
band of spatial frequencies u, whereas instabilities of other explicit schemes occur, typ- 
ically, for all frequencies higher than a certain minimum value. Let us note a curious 
consequence of this spectral selectivity of the instability of the split-step method: The 
method can be found stable even when Az exceeds the threshold (11. 7p . This can occur 
when the width of the instability peak is less than the mesh size Ao; in the frequency 
domain [6]; see also Appendix 3. 

In this paper, we will show, first via numerical simulations and then by an analyt- 
ical calculation, that the properties of the instability of the split-step Fourier method 
on the background of a soliton solution are considerably different from properties of 
such instability on the monochromatic wave background. In particular, these new 
properties are even more distinct from some propeties of instabilities of well-known 
finite-difference explicit schemes. Highlights of these new properties include: high 
sensitivity to (i) the step size Az and (ii) the length T of the time window (recall 
that time t in (11.11) is a spatial coordinate), and also (iii) a violation of the so called 
principle of "frozen coefficients" (see below). Property (i) has been considered both 
numerically and analytically for the monochromatic wave background, and numeri- 
cally for the soliton background, in [8]. We will explain later on that this property, 
i.e. the high sensitivity of the instability to the step size Az, has different dependence 
on the time window length T and also requires a different mathematical description, 
for the monochromatic and soliton backgrounds. As for properties (ii) and (iii), they, 
to the best of our knowledge, have not been previously systematically studied for any 
numerical scheme. 

In Section 2, we will present results of our numerical simulations of Eq. (11.11) by 
the split-step Fourier method (II. 2p on the background of the soliton 

u,oi = A y^sech (^^) exp [iKz] = U{t)e'^\ K = A\ (1.8) 
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These will illustrate the unusual instability properties listed in the previous paragraph. 
In Section 3, we will develop an analytical theory of the instability of the slit-step 
Fourier method on the background of the soliton which explains all these observed 
propeties. A comparison of this theory with numerical simulations is presented in 
Section 4. In Section 5 we summarize this work. Appendices 1 and 2 contain auxiliary 
results, and Appendix 3 recovers the results of Weideman and Herbst [6], illustrated 
in Fig. [H via the analytical method presented in Section 3. 



2 Numerical study of the instability of the split- 
step Fourier method on the background of a soli- 
ton 

We numerically simulated Eq. (11. ip with /3 = — l,7 = 2on the background of the 
soliton (II. 8p with A = 1 using algorithm (11. 2p . We added to the initial soliton a very 
small white (in t) noise component, which served to reveal unstable Fourier modes 
sooner than if they had developed from the round-off error. Thus, the initial condition 
for our numerical experiments was 

Uo{t) = A sech (At) +^{t), A = l, (2.1) 

and ^(t) was a Gaussian random process with zero mean and the standard deviation 
10-10. 

To begin, we considered the spatial grid — T/2 < t < T/2 with 2^^ grid points 
and width T = 327r, i.e. about two orders of magnitude wider than the soliton. This 
results in the spectral grid being the interval —32 < u < 32 with the frequency 
spacing of Au = 2tt /T = 0.0625. As expected, we did not observe any instability for 
Az < Azthresh ~ 0.0031. Abovc the threshold, we ran the simulations with Az ranging 
from 0.004 to 0.006 with the increment of 0.0001 (i.e., Az = 0.0040, 0.0041, . . . 0.0060) 
up to the maximum distance of Zmax = 500 and observed instability only at a few values 
of Az listed in Table 1. A typical Fourier spectrum of the unstable solution at z^^^^ is 
illustrated in Fig. |2]for Az = 0.0040. The instability increment (see ( ll.Sbl) ) listed in 
Table 1 was computed as 

(peaks' exponent) — (noise's exponent) 
Re(A) = In 10 , (2-2) 

where the peaks's exponent refers to the average of the decimal exponents of the two 
peaks in Fig. Et^b), and the exponent corresponding to the average noise level was 
estimated to be —8.8 for this set of simulations. Estimate (12. 2p may not be very 
accurate, as the so computed increment depends on the noise levels at the frequencies 
of the unstable Fourier modes (hence the peaks in Fig. [2] are slightly different), but it 
still provides a reasonable measure of the instability rate. 
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Az 


icrt 




rignt 


Re(A) 

V / 


0.0040 


-0.72 


28.03 


0.66 


0.019 


0.0048 


-0.39 


25.58 


0.35 


0.031 


0.0054 


-0.62 


24.12 


0.57 


0.022 


0.0055 


-0.52 


23.90 


0.48 


0.024 


0.0058 


-0.46 


23.27 


0.42 


0.022 



Table 1: Parameters of the unstable frequencies' peaks when T = 327r, number of grid 

points is 2^°, and Az is varied as 0.0040, 0.0041, 0.0042, . . . ,0.0060. The notations 

,(+) 

right, left 



'^right left ^'^^ introduced in Fig. [2l 



From Fig. |2] we first observe that the instability peaks at negative frequencies are 
not refiectionally symmetric relative to such peaks at positive frequencies, in contrast 
to the instability peaks on the background of a monochromatic wave (see the caption 
to Fig. [1]). Rather, for all the simulations we ran, the negative-frequency peaks appear 
to be a shifted replica of the positive-frequency peaks. Moreover, the frequencies of the 
left and right peaks on the same side of u = are slightly asymmetric with respect 
to Uj^: see Table 1. To further investigate how the peaks are related to one another, 
we ran a simulation where we placed a narrow filter that totally suppressed the field 
around uj = oo^^^^^- As a result, the peak at a; = 0;^;"^^ was suppressed but the peaks at 
00 = ool^^ remained intact. This corraborated the aforementioned observation that the 
pairs of peaks appear to be shifted replicas of each other rather than refiected about 
oj = 0. In light of this, it might be somewhat surprising that the peak's frequencies are 
found to be related by a refiectional symmetry: 

^left - ~^right' ^right " ~^left • 

All these observations are explained by the theory in Section 3. 

Another conspicuous observation, made from Table 1, is that there seems to be 
no apparent order in the frequencies and heights of the instability peaks as func- 
tions of Az. To expand on this, we zoomed in on the interval 0.0045 < Az < 
0.0046, at the end points of which we had found no instability. We varied Az = 
0.00451, 0.00452, . . . 0.00459 and found that the instability occured at Az = 0.00451 
and persisted up to Az = 0.00456, with its increment Re(A) gradually decreasing from 
0.036 to 0.014 and the inter-peak spacing u}^i~gu~^Mt gradually increasing from 0.25 to 
1.56. A weaker instability, with the inter-peak spacing of about 2.0, was also observed 
at Az = 0.00459; note that no instabihty was observed at Az = 0.00457 and 0.00458. 

One could argue that such an irregular dependence of the instability on Az occurs 
simply because the width of the band of unstable modes is just slightly greater than 
the mesh size Au of the frequency grid. Indeed, the unstable band's width estimated 
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Figure 2: The spectrum of the solution of ([EI]) and at z = 500 with Az = 0.0040, 
computed by method (ll.2p . Other parameters are hsted in the text after Eq. (12. ip . 

from the monochromatic-background case is 7A^/(|/3|(X',r) ~ 0.07 for the parameters 
hsted above (see the caption to Fig. [1] and formula ( 1A3.6I) in Appendix 3), while 
Au = 271 /T = 0.0625. In such a case, the instability features should be strongly 
sensitive to where inside the instability band the frequency 2iTi/T falls. This was 
pointed out in [8]; see also the end of Appendix 3 below. Then one would expect 
that the aforementioned sensitivity is to be alleviated by taking a wider time window, 
because then the frequency mesh size Au would decrease and more frequencies from 
the numerical grid would fall into the unstable band. 

To verify the validity of such an argument, we re-run the simulations described 
above taking a four times wider time window, i.e. T = 1287r, while retaining the same 
At (thus quadrupling the number of grid points). For these new simulations, we esti- 
mate that instability peaks should contain about four grid points and hence may expect 
that the high sensitivity to the value of Az reported above is to be alleviated. What 
we find instead is an opposite of this statement: the dependence of the locations and 
heights of the instability peaks remains at least as irregular as for the smaller value of 
T, but now the instability is observed "more often" than in Table 1. In Figs. |3] and H] we 
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Figure 3: Instability increment and unstable peaks' half-separation at z = 2000, as a 
function of the step size Az. The initial condition is (12. ip . The numerical domain has 
length T = 1287r and contains 2^^ grid points. Note that for Az = 0.0040, where no 
instability is found, we defined the corresponding Unght — <^\cft = 0. The open circles 
and solid line correspond, respectively, to the numerical results and to the results of 
analytical calculations reported in Section 4.2. 

plot the observed values of the instability increment and half of the inter-peak spacing, 
(^Sht - ^icft)/2, when Az is varied between 0.0039 and 0.0050 with step 0.0001 and 
between 0.00471 and 0.00479 with step 0.00001. Since the instability rates have now 
been found to be about a factor of four lower than in Table 1, we had to run our simu- 
lations up to a greater distance, Zmax = 2000. Note that the instability characteristics 
reported for 0.0040 < Az < 0.0050 in Table 1 do not match those shown in Fig. [31 
Moreover, we find that the spectra of unstable modes may look qualitatively different 
than in Fig. d Namely, these spectra for Az = 0.0044, 0.0045, 0.00474, 0.0049, look 
like the one shown in Fig. El^a), while for Az = 0.0050 it is shown in Fig. El^b). Also, 
contrary to our expectation, we observe that in most cases the instability peaks still 
contain only one grid point; exceptions are the central peak for Az = 0.0050 and the 
peaks for Az = 0.00478, which are spaced very closely. This fact will be emphasized 
when we describe a challenge in the analytical description of the instability in Section 
3. 

We have already mentioned that keeping the spatial mesh size At intact but increas- 
ing the spatial window's length T by a factor of four considerably changes parameters of 
the instability. We now show that just slightly changing T (and hence correspondingly 
slightly changing At) may change instability parameters dramatically. Such a sensi- 
tivity to the length of the spatial window is not observed for finite-difference methods. 
Note that T = 1287? ^ 402.1. In Table 2 we list the instability parameters observed 
when we decrease T by about 1% or even less. Again, these results appear to be totally 
irregular. Similarly to the situation mentioned at the end of the previous paragraph. 
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Figure 5: "Less typical" instability spectra. Both panels are for initial condi- 
tion fl2.ip . the propagation distance z = 2000, and the numerical domain with 
T = 1287r and 2^^ grid points. (a) Az = 0.0049; similar spectra are found for 
Az = 0.0044, 0.0045, 0.00474; (b) Az = 0.0050. 



only the closely spaced instability peaks at T = 398 contain several grid points; the 
peaks for the other values of T each contain only one grid point. 

Finally, let us show that the so called principle of "frozen coefficients" [9], which 
is known to apply to finite-difference schemes, does not hold for the split-step Fourier 
method (11. 2p on the background of a soliton. This principle says the following. Sup- 
pose one has an evolution equation with a spatially varying coefficient, say c{t) (recall 
that in this paper, t is a spatial variable, while z is the evolution variable). Near each 
t = to, this coefficient can be approximated by c{to). Then one can apply the standard 
von Neumann stability analysis to the equation with the constant ( "frozen" ) coefficient 
c(to)- Then the scheme is deemed unstable if such an analysis reveals instability for at 
least one value of c(to)- This principle works quite well for finite-difference schemes. 
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T 


U/i^ff — 0/77- 


rigiit " 


Re(A) 


397 


-0.44 


0.40 


0.0062 


398 


-0.10 


0.08 


0.0096 


399 


-1.09 


1.02 


0.0024 


400 


-0.33 


0.29 


0.0076 


1287r 


-0.39 


0.36 


0.0076 



Table 2: Parameters of the unstable frequencies' peaks when /S.z = 0.0043 (cUyr = 27.03), 
number of grid points is 2^^, and T is varied. 

An intuitive explanation fot this is that the unstable modes usually have high spa- 
tial frequency, and so they "see" a relatively slowly-varying coefficient c(t) as being 
approximately constant near each t = to. For the split-step Fourier method, the un- 
stable modes also have high spatial frequencies, ~ ±00.,^. Yet, the principle of "frozen 
coefficients" does not apply to this method, as we illustrate below. 

Indeed, we have already seen that the instability of method (11. 2p on the background 
of a soliton is different from that on the background of a monochromatic wave of the 
same amplitude. This fact was originally stated by Weideman and Herbst (see the 
last section in [6]). Now we will present three examples that deal with a multi-soliton 
background. Conclusions from these examples are at odds with our intuition based 
on the experience with finite-difference schemes, and they reveal yet another way in 
which the principle of "frozen coefficients" is violated for method ( 11. 2p . As at the 
beginning of this Section, let us use the numerical domain with T = 32-7? ^ 100.5 and 
210 gj^jj^ points, but instead of the sm^^/e-soliton initial condition (12. ip . consider three 
well- separated solitons: 

uoit) = sech {t + 33.5) + sech (t) + sech {t - 33.5) + ^{t), (2.4) 

where ^(t) is the same as in (12. ip . The spacing between the solitons is chosen to be 
sufficiently large so as to avoid their interaction. According to the principle of "frozen 
coefficients", the instability with this initial condition must be the same as with (12. ip . 
However, it turns out to be quite different. First, for Az = 0.0040, when there was an 
instability (see Table 1) for the initial condition (12. ip . there is no instability for initial 
condition ( 12. 4p . Second, for Az = 0.0052, when there was no instability according to 
Table 1, now there is a strong instability with the increment of about 0.10 (i.e., more 
than three times greater than the largest increment in Table 1). Third, for Az = 0.0058, 
the instability was observed with both initial conditions ( 12. ip and (12. 4p . but for the 
latter case it was more than three times as strong (having the increment of 0.072). 

To summarize, in this Section we have presented results of numerical simulations 
of the split-step Fourier method (11. 2p with soliton initial conditions. These results 
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illustrate the following features that are drastically distinct from features exhibited by 
unstable finite-difference schemes: 



• The spectra of the unstable modes exhibit strong and irregular sensitivity to the 
values of the step size /\z and the length of the spatial window T; 

• Instability on the background of several well-separated — and hence non-interacting 
— pulses can be drastically different from the instability on the background of a 
single pulse (or a different number of pulses). 

The theory that we will present in the next section is able to quantitatively explain all 
these features. 



3 Analytical theory of the instability of the split- 
step Fourier method on the background of a soli- 
ton 

3.1 Key idea and main challenge 

First, we will show that a straightforward application of the von Neumann analysis 
is unlikely to provide analytical insight about stability or instability of the split-step 
method on the soliton background. Then we will outline the idea of an alternative 
approach and point out the mathematical challenge that such an aproach would have 
to resolve. 

Substituting expression f ll.Sap with u = Usoi into algorithm (11. 2p and keeping only 
terms linear in «„, one obtains after taking the Fourier transform: 

J-^+i] =e^^-'^-^[e*^l"-l'^^(n„ + ^7A^(n2„i< + |«soi^^ , (3.1) 

where Ugoi is given by (II. 8p . The right-hand side of (13. ip describes couphng of Fourier 
modes J^[un]{uji) with different ui = 2TTi/T, i = 0, ±1, ±2, . . . via, e.g., the convolu- 
tion term J^^I-UsoipMn] , because |msoi| is a function of the spatial variable t. (Note that 
this problem did not occur in the stability analysis on the background of a monochro- 
matic wave (II. 4p since there l^l = ImcwI is t- independent and hence the corresponding 
equation, studied by Weideman and Herbst, coupled only the modes ui and —Ui.) 
Solving for eigenvalues of such a coupled multi-mode system, while numerically fea- 
sible, would unlikely provide any insight of how the instability can occur. Such an 
insight is provided by an alternative analytical approach described below. 

The key idea of this approach comes from the instability spectra shown in Figs. [H |2l 
and|5l Namely, we note that only narrow bands of Fourier modes with frequencies near 
iu^TT can become unstable. The reason for this will be presented later. Therefore, it is 
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appropriate to seek the numerical error Un, defined in (11 .Sap , as consisting of two quasi- 
monocfiromatic waves whose carrier frequencies are approximately iw^. These waves 
can become unstable via their interaction mediated by their scattering on the soliton. 
Since these are high-frequency waves (see (ll.6p and note that Az ^1), they "see" the 
soliton, whose temporal width is 0(1), as a narrow, and hence small, perturbation. 
This is one of the reasons that the instability is weak, as seen in Section 2. Using the 
weakness of the instability, we will approximate its evolution by a differential, rather 
than difference, equation. Analysis of the former kind of equation is considerably easier 
than that of the latter one. 

The main mathematical challenge that this approach needs to address is a three- 
scale nature of this problem, as can be seen from Fig. |2l One scale is set by ^ 1, 
where the strong inequality is a consequence of Az <^ 1 (see (II. 6p ). as is usually the 
case in simulations. Accordingly, let us define a small parameter 

e= — <1, e = 0(VA^), (3.2) 

which will play a prominent role in what follows. In addition to this 0(e^^) scale, 
there is also a scale 0(1): The separation between u-,, and the frequencies of the most 
unstable modes, seen as peaks in Fig. |21 appears to be of this order of magnitude. 
Finally, the third scale is set by the spectral width of the instability peaks. Indeed, as 
noted in Section 2, even for the highest spectral resolution reported there, those peaks, 
in most cases, contain only one grid point. (The pedestals seen around these peaks can 
be shown to occur due to modulational instability on the background of the peaks, i.e. 
are not directly caused by a numerical instability.) Based on the information provided 
by Figs. [2] and O it is even impossible to tell whether this third scale is determined by e. 
We will show below that it is rather determined by several parameters of the problem. 
Incidentally, let us note that the instability on the background of a monochromatic 
wave has only two scales: the location ±Utj = 0(e~^) of the instability peaks and the 
peaks' width, 0(e); see Fig. [1] and Appendix 3. 

3.2 Details of the theory 

First, let us simplify the right-hand side of (13. ip by noting that in practice, Az is always 
taken so as to guarantee 'y\usoi\'^Az <^ 1. Then, discarding terms 0[{'-y\usoi\'^Az)^) , one 
reduces (13. ip to 

J'iun+i] = e*''"'^' ^ [un + i-fAz{uliul + 2\u,ol\^Un)] • (3.3) 

Let us now use the observation made in the previous subsection: As follows from Figs. 
[H 121 andini the instability occurs only in narrow spectral bands near icUvr, ±uj2tt, etc. 
We will focus on the instability near ±Ut^] the analysis near ±U2tt etc. is similar. (In 
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Appendix 2 we will show that the instability can occur only near ±u;,r, ^^2n, etc., 
thereby recovering the results presented below.) Accordingly, let us rewrite (13. 3p as 

= (-1) ■ e^/3(-'-^)^^ ^ [m„ + tjAziulK + 2|M,oipM„)] . (3.4) 

where the (—1) in front of the right-hand side occurs due to the definition of u^^, 
Eq. dLS]). Note also that (3{u'^-ul)Az = 0(e), which follows from ([32]) and the fact 
that we consider frequencies satisfying — {±ijj.j,)\ = 0(1); see the end of Section 3.1. 
Thus, exp[i/3(u;^ — u^Az] = 1 + 0(e) ~ 1 for the frequencies of interest. 

To enable further analysis of the still intractable difference equation (13.41) in the 
frequency domain, let us convert it into a partial differential equation in the time 
domain. To this end, we first observe that for the variable 

C„ = (-1)"m„, (3.5) 

Eq. ( 13. 4p describes a small increment occurring from nth to {n + l)th step: 

^[^„+i] = e*/3(-'-^)A- jr + t^Aziul,v: + 2\u,oi\''vn)] ■ (3.6) 

Indeed, from (13.61) and the estimate exp[i/3(u;^ — uDAz] = 1 + 0(e) it follows that 
Vn+i — Vn = 0(e). Accordingly, we define a continuous variable v{z,t) which at Zn = 
nAz is related to {t„ by (13. 5p . Then, we show in Appendix 1 that a variable w = v e~^^^ 
satisfies an equation 

= -i(3{wtt + ujIw) - iKw + i-fU^{w* + 2w), (3.7) 

where K and U = U{t) are defined in (II. 8p . and we have neglected terms of magnitude 
0(e) (see Appendix 1). Note that the first group of terms on the right-hand side of 
(13. 7p has the order of magnitude 0(l/e) (see the text after (13.40 ). and thus it describes 
waves rapidly oscillating in both z and t, as we have announced in Section 3.1. 
We seek a solution of (13. 7p as 

w = pe-*^^* + m*e^'"^*, (3.8) 

where p{t,z) and m{t,z) are assumed to vary in time on a scale 0(1). This agrees 
with the numerical observation in Section 2 that frequencies of unstable Fourier modes 
differ from ±u;^ by 0(1) in most cases. Substituting (13. Sp into (13. 7p and separating 
terms proportional to exp(±ici;^t), one obtains: 

= -2{f3/e)pt - iKp - if3pu + t^U\m + 2p), (3.9a) 

777-2 = 2((3/e)mt + iKm + il3mtt — i'jU'^{p + 2m), (3.9b) 

where we have used definition (13.20 . Using again the aforementioned numerical obser- 
vation from Section 2, we seek solutions of (13. 9 p in the form 

P = PsIow(t) exp [-int + 2i{f3/e)nz + f3Az] , (3.10a) 
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m = msiow(T) exp [iQt + 2i{/3/e)Qz + /3Az] , (3.10b) 

where 

n = 0{l) (3.11) 
labels a particular Fourier mode of w, parameter A is to be determined later, and 

T = et, A = Ar + 'iAi. (3.12) 

Note that /3A/j equals Re(A) in ( ll.Sbl) . When writing that psiow and msiow, as well as 
P(r) and M(r) in f l3.15p below, are functions of the "slow" time r, we mean that 

(Psiow)t, ("^siow)t, Pt, Mt are all of order 0(e). (3.13) 

The order-of-magnitude estimate (13.111) . which we empirically deduced from numerical 
experiments, will be generalized at the end of Section 4. 

Substitution of (13.101) into (13. 9p produces a pair of ^-independent equations: 

2(1 + ef^)(psiow)r = {-^K//3 + t^^ - A) p.iow + t{i / I3)U\t / e) (msiow e'^^^/^ + 2j9,iow) , 

(3.14a) 

2(l-e^])(mslow)r = {-iK/P + i^"" + A) m,,,^+i{i / f3)U\r / 1) (psiow e"'*^^/^ + 2m,io„) . 

(3.14b) 

In writing these equations we have neglected terms O(e^), which will be justified by 
subsequent calculations. Note also that the soliton background, f/^(r/e), presents a 
narrow obstacle for PsIow(t) and msio^ir). Since outside the soliton, Ps\ow and rrisiow are 
not coupled to each other, we use yet another substitution: 



Psiow = P exp 



"^siow = M exp 



i(-ir//3 + _ A 

2(r+d7) ^ 
'i{-K/P + ^l'')+A 



T 



2{l-eVL) 

Then P and M change only in the vicinity of the soliton according to 



(3.15a) 
(3.15b) 



= ^^l^'l'"/^)^ ^2^n■r/e+Oil)■r ^ , (3.16a) 

= '^^^'^^/'l 2^n■r/e+Oil)■r gM) . (3.16b) 

2/3(1 -en) ^ ' ^ ^ 

Integrating these equations over the entire time window — eT/2 < r < eT/2, one 
obtains: 

P(+eT/2) - P(-eT/2) = ^ {J^[U^]{-2n) M(0) + 2J^[U^]{0) P(0)) , (3.17a) 

M{+eT/2) - M(-eT/2) = «e (J^[t/2](2fi) P(0) + 2T[U%0) M(0)) . (3.17b) 

Zp 
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Here we have neglected terms O(e^) and used the definition ( II. 3p of the Fourier trans- 
form and the fact that the sohton is centered at r = 0. These relations prompt us to 
denote relative jumps eJp and eJu that P and M undergo across the soliton: 

P(±eT/2) = P(0) (1 ± eJp/2), M(±eT/2) = M(0) (1 ± eJu/'i). (3.18) 

Thus, by definition, Jp^u = 0(1). According to f l3.17p and with the same accuracy, 
these jumps satisfy 

Jp = i^ {j^[U^]{-2n) R + 2J^[f/2](0)) , (3.19a) 

7 ^^^.,2v.^^ 1 



JM = ^^[J'mi2n)- + 2T[U']iO)j, (3.19b) 

where R = M(0)/P(0). 

Let us pause for a moment and recall what we are trying to do: We want to deter- 
mine the instability increment /3Ar for the Fourier mode whose frequency is related to 
Q by fl3.10p and (13.80 . The two equations ( I3.19P for three unknowns Jp,M and R are 
insufficient for this purpose; note that they do not even involve A. The missing rela- 
tions that will allow us to complete our task are supplied by the periodicity condition 
satisfied by P{t) and M(r) and are obtained as follows. First, note that the numerical 
error Unit), and hence w{t,z) defined before (13. 6p . satisfies periodic boundary condi- 
tions in t by virtue of the split-step method (11. 2p using the discrete Fourier transform. 
Second, since p(t,z) and m*{t,z) have different ^-dependences (see (I3.10p ). then by 
virtue of (13. 8p each one of these functions must satisfy periodic boundary conditions 
in t. Third, along with I^M>, dSJO]), and f l335|) . this implies that 



P(+6T/2) _ 
P(-eT/2) " """"P 

M(+eT/2) 
M(-eT/2) 



exp 



^(rTrf]^^^+^^"-+''^^. 

el + i[uJt, — 



(3.20a) 



(3.20b) 



2(1 -efi) 

Finally, using the identity exp[2i7rA^] = 1 for integer and Eqs. (13.200 and (I3.18p . 
and neglecting terms O(e^), we obtain: 

[i {K/13 - n^) + A] (1 - en) eT/2 + i{6u^ + n)T = 2mNp + eJp, (3.21a) 

[i {K/P - fi2) - A] (1 + en) eT/2 + i{6uj^ - fi)T = 2z7rA^M + e^M, (3.21b) 

where Np^M are some integer numbers and Su-j^ is "the fractional part" of w^: If 
= 27m /T (where n is not necessarily an integer) and N.,^ is the integer part of n, 
then 6co^ = 2Tx{n - N^)/T. Note that 5u^ = 0{1/T). 

We can now justify neglecting the terms O(e^) in (I3.14p and in subsequent calcula- 
tions. Indeed, if such terms had been retained, they would have contributed amounts 
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O(e^) and 0{e^Q ■ eT) to fl3.2ip . The former amount would be a higher-order contri- 
bution than that provided by terms eJp^M, which we need to determine. The latter 
amount, strictly speaking, depends on the order of magnitude of (eT), but it will be 
clear from our subsequent calculations that even terms 0{eQ ■ eT) can be neglected in 
the leading-order analysis. 

We will now use Eqs. (13.191) and (I3.2ip to determine for what values of Q the 
instability increment (SAn is nonzero. To this end, we first subtract Eqs. (13.2ip from 
each other and take the real part, obtaining: 

Next, adding Eqs. (I3.2ip one obtains: 

e{Jp + Jm) = {iK/(3 - iVl'^ - A en) eT + 2i5uj^T - 2i7i{Np + Nm) ■ (3.23) 

Using Eq. (I3.22p one notices that the real part of the the right-hand side of (I3.23P is of 
order O(e^) and hence should be neglected. Thus we conclude that in the main order, 
( Jp + Jm) is purely imaginary. For future use, we also display the result of taking the 
imaginary part of the difference of the two equations (I3.2ip : 

Im( Jp - Jm) = (- {K//3 - n^) eVl + Aj) eT + 2VLT - 2-n{Np - Nm) ■ (3.24) 

Since Ap is proportional to the real part of (Jp — Jm), we solve for the latter 
quantity using Eqs. (I3.19p . To that end, we first solve for R by adding these equations 
and then substitute the answer in their difference, obtaining: 



Jp - Jm = ±j\J (2^[f/2](0) + z(/3/7)(Jp + Jm))' - |-F[t/2](2fi)|' . (3.25) 



Now recall that, as noted after f l3.23p . (Jp + Jm) is purely imaginary. Also, J-'[f/^](0) 
is real. Then the real part of the right-hand side of ( 13.25P is nonzero when 

- |J^[f/'](2^])| < 2J^[f/2](0)+i(/3/7)(Jp + Jm) < \J'[U\29)\ . (3.26) 

Under this condition, one also has 

Im( Jp - Jm) = 0. (3.27) 

Thus, the instability increment, /3Ap, is found from Eqs. (13.221) and fl3.25p . where 
(Jp + Jm) is determined from (13.230 . The last two equations, in their turn, involve 
three yet undetermined quantities: f2, which labels the frequency of an unstable Fourier 
mode, parameter Aj introduced in (I3.10p . and (A''p + A^m)- We will now show that 
within the accuracy adopted in our calculations, the former two quantities enter all 
equations in a unique combination {VL + eA//2), and hence the number of yet unde- 
termined quantities reduces to two. Indeed, upon substitution of (13.150 into (I3.10p . 
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we observe that up to terms O(e^), both t- and 2;-dependences of p and m involve Q 
and A/ only in the aforementioned combination {Q + eAj/2). Next, by inspection of 
formulae (13.231) and f l3.24p . one can easily see that, within the same accuracy, they also 
involve Q and A/ only in that combination. This means that one can set A/ = 0. Then 
from (13:21) and (13:271) one finds: 

n = 7r{Np-NM)/T, (3.28) 

where we have discarded terms O(e^). Then (13.231) is rewritten as 

^(Jr + Jm) ^ - (f + T + (fil^L^ ^ 2MN. + Nm) ^ ^^^^^ 

Finally, one substitutes the last two equations into (I3.26P and determines those values 
of {Np ± Nm) where the instability can occur. The instability increment is computed 
from ( I3.22P and ( I3.25p . and the frequency of the unstable mode, from ( 13. 28 p . Examples 
of such a calculation, producing the theoretical results shown in Figs. [3] and HI are given 
in Section 4.2. 

In Appendices 2 and 3 we present related technical results. Namely, in Appendix 2 
we show that if one seeks unstable modes not specifically near ^Ut^, as we did at the 
beginning of this subsection, but instead near an arbitrary pair of frequencies iwo, one 
discovers that the instability can arise only near ±uJt,, ±uj2n, etc. In Appendix 3 we 
modify the analysis of this subsection to apply to a monochromatic-wave background 
( II. 4p and thereby recover results that can be deduced from those obtained by Weideman 
and Herbst ^6j. 

4 Validation of the theory 

In the first subsection below, we will give qualitative explanations of the instability 
features described in Section 2: the locations and widths of the instability peaks, and 
high sensitivity of the instability to the step size, the time window length T, and 
the shape of the background solution. In the second subsection, we will work out 
an example showing how Eqs. (I3:22|) . ([S^SD, ( K26\f . (I3:29|) . and ^28\f were used to 
compute the increment and frequency of the unstable Fourier modes reported in Figs. [3] 
and m We will conclude by generalizing the order-of- magnitude estimate ( 13. lip for the 
unstable modes' separation frequency, 2fl. 

4.1 Qualitative explanation of instability features reported in 
Section 2 

The results of Section 3.2 allow us to explain why the instability peaks, shown in Fig. 
|2l are not reflectionally symmetric about u = (see the paragraph after Eq. ( 12. 2p ). 
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Indeed, from Eqs. (13.81) . ( I3.10p . and (IS.lSp . one sees that the frequencies of two coupled 
unstable modes are 

{uj^ + n- e{-K/l3 + n^)/2) and - {uj^ - Q - e{-K/l3 + n^)/2) . (4.1) 

Thus, given the sign difference of the second terms inside the parentheses above, the 
mode at uj^^^^^_ is coupled to the mode at weight ^^'^ '^^^ tliaX at oj[^^^, as it would have 
been in the case of refiectional symmetry. 

Similarly, two other features of the instability spectra reported in Section 2 can 
also be explained. First, note from (I3.28P and (I3.29p that if a value > is found 
to correspond to an instability peak, then so is —VL < 0. This observation, along with 
relations 

4St = ±{[io^~e{-K/P+n^)/2]±n) and 4tt = ±{[io^-e{-K/P+n^)/2]Tn), 

(4.2) 

which follow from (14. ip . explains why relations (12. 3p hold. Second, the slight asymme- 
try of the frequencies weight u^^^ about the respective ±0;,^, is also easily explained. 
Indeed, the frequencies of, say, the peaks at w^j^ht ^icfh ^^^^^ from (14. 2p . are 
centered about [wtt — e{~K/f3 + i7^)/2] and not about u,^. The maginitude of the 
shift, —£{—K/l3 + Q'^)/2, agrees with the numerically observed values. For example, in 
the experiments reported in Table 1, K = 1, P = —1, e ~ 0.04, and Q ~ 0.5, and hence 
—e{—K/ (3 + Q?)/2 ~ —0.025. This should be compared to the experimental values of 
[('^left ~ '^1-) + ('^right ~ '^1')] /2) which from Table 1 are seen to vary between —0.02 and 
-0.03 0. 

Next, we can explain why in most mentioned in Sections 2 and 3.1, the 

instabihty peaks contain just one nodjfl. Consider Eqs. (I3.28P and (I3.29P and assume 
that at frequency labeled by VLq, corresponding to a particular {Np — Nm)q, there is 
an instability. The frequency at the adjacent node differs from this Qq by 27r/T, and 
hence, according to (13.280 . the corresponding (Np — Nm) differs from {Np — Nm)o by 2. 
Consistently with this, one can take (Np+Nm) = {Np+Nm)o- Then the corresponding 
value of the right-hand side of (I3.29P differs from that value at Qq by 2n^{Np—NM)o/T ■ 
2 = AtcQq. On the other hand, the interval of values where i[Jp + J^) corresponds 
to an instability is found from (I3.26P to have the width 2 |(7//3)J^[?7^](2fio)|- Thus, 
whenever 

Ann > 2\{^//3)J^[U^]{2n)\, (4.3) 

more detailed comparison would require keeping at least one more significant digit in the data 
of Table 1, but such a comparison does not appear to be needed. Rather, it is the quantitative 
agreement of our theory and numerical experiments, reported in Figs. [3] and HI and presented in the 
next subsection, which seems to be the most important test confirming the validity of the theory. 

^As we already mentioned in Section 2, the pedestals around the peaks arise due to a non-numerical 
— modulational — instability, and those pedestals are hence unrelated to the foregoing explanation. 
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the instability peak can contain only one node. Using expression (II. 8p for U{t), one 
finds that condition (14. 3 p holds for Q > 0.46. In other words, it is only when the peaks 
at cjieft and Wright are separated by less than 0.92 that they can contain more than one 
node. In practice, the separation should be even smaller, given a random location of 
the node with the unstable frequency within the instability interval. In fact, among 
the simulations reported in Section 2, we observed multiple nodes per peak only when 
(wright — i^ieft) was about 0.25 or less. 

The high sensitivity of the instability to the length T of the time window, which was 
highlighted in Section 2 (see also Table 2 there), is also easily explained using (I3.29p . 
Suppose this length is changed from Tq to T, so that the relative change (T — Tq)/Tq 
is small. Then the right-hand side of (I3.29P is changed by an amount 



The coefficient in front of (T — Tq)/Tq is large due to Tq being large and the terms in 
parentheses being 0(1). Note that this does not necessarily imply that the instability 
will disappear, since new values of Np and Nm may exist for which the right-hand 
side of (I3.29P will be inside the instability interval (I3.26P : see the next subsection for a 
quantitative example in a similar situation where not T but Az, is varied. 

Finally, while instability on the background of several well-separated solitons cannot 
be computed from Eqs. (I3.25P and (I3.26p . which were obtained for a single soliton, one 
can still explain why the instability is sensitive to the number of the solitons and their 
relative locations. Indeed, the number of solitons determines the number of jumps in 
the "slow" variables P and M; see (I3.16p - fl3.19p . Its is those jumps that allow modes 
that exponentially grow/decay (in t) away from the jumps to exist in the presence of 
periodic boundary conditions in t; see fl3.2ip . And its is those modes, exponentially 
growing or decaying (as opposed to purely oscillating) in t, that are also exponentially 
growing (i.e., are unstable), in the evolution variable z; note the same exponent A 
in (I3.15P and in fl3.10p . As for the relative locations of the solitons, they (as well as 
the solitons' phases) determine the counterparts of parameter R in a generalization of 
Eqs. f l3.19p for a multi-soliton background. 

4.2 Example of calculation of instability increment and fre- 



Here we will explain in detail how the instability increment and frequency for the 
generic case are found. Such a generic case corresponds to, e.g., Az = 0.0043 in Fig. |3l 
for which we will present the calculations below. Its instability spectrum is qualitatively 
the same as that shown in Fig. O Then we will comment on less common cases, whose 
spectra are shown in Fig. [51 We will conclude this section by establishing a general 
dependence of clS cL function of e and T. 




(4.4) 



quency 
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For all cases considered here, the parameters in the background soliton (ll.Sp are: 
P = —1, 7 = 2, A = 1, whence K = 1 and U{t) = sech(t). The time window is 
T = 1287r and the number of grid points is 2^^. These are the same parameters as were 
used in obtaining Figs. |3]-|5l 

To begin, let us recall from f l3.25p and f l3.26p that for the mode labeled by Q to be 
unstable, the left-hand side of fl3.29p must fall within the interval 



(4^[f/2](0) - 2|^[f/2](2fi)| , 4^[f/2](0) + 2|^[f/2](2f])| ) 



(4.5) 



(here we have used that 7//? = —2). While f2 is yet undetermined and hence |^[f/^] (2f2) | 
is not known exactly, it cannot exceed J-'[f/^](0) = 2, and hence interval (14. 5 p is always 
inside the interval (4, 12). Thus, this interval contains values of order 1, i.e. much 
smaller than any of the terms on the right-hand side of (I3.29p . which are proportional 
to the large parameters T and 1/e = u.,,. This observation motivates our strategy of 
finding suitable values of {Np ± Nm)- 

First, given the above values of parameters, one finds u-,, = 1/e ~ 27.03, 6ut, ~ 
0.0140, and then the first term on the right-hand side of (I3.29p : 



K 26u^ 
\ 

(3 e 



T ^ 96.86. 



Next, since the first term in the second parentheses, {n^Np — Nm)Y /T, 
positive, we seek such an integer value for {Np + Nm) that an expression 



K 25uj^ 
\ 

/3 e 



T + 



27i{Np + Nm) 



(4.6) 
is always 

(4.7) 



is as close to zero as possible but negativ^. Since 27r/e ~ 169.84, the corresponding 
Np + Nm = —1, yielding the value of —72.98 for (14. 7p . Finally, an integer value of 
{Np — Nm) is sought to make the value of the right-hand side of (I3.26P fit within the 
interval (14. 5p . To this end, one first estimates {Np — Nm) by forcing that right-hand 
side to equal zero: 



(Np - Nm) 



M I estimated 



K 26uj^ 
\ 

(3 e 



T 



271 {Np + N, 



M) 



(4. 



Here (A''p — A^jv/) estimated is not an integer; e.g., for the example above it is approximately 
54.53. Then a nearby integer value for {Np — Nm) is found by trial and error so that it 
makes the right-hand side of (13.291) fit within the interval (14. 5p . For the above example, 
this value is {Np - Nm) = 57, yielding Q ^ 0.44 (see ([32HD) and PAp = 0.0067 (see 
(I3.25P and (I3.22p ). Comparison of so computed f3Ap and Q with respective values 
found numerically is shown in Figs. [3] and |4] for various values of Az. 



^Modifications of tliis requirement will give rise to "less common" cases, considered later. 
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Two notes are in order about the last step described above. First, since both Np 
and Nm are integers, then {Np + Nm) and {Np — Nm) are either both odd or both 
even. Second, once a guess is made about {Np — Nm), then Q is computed from fl3.28p 
and hence |j^[f/2](2fi) | in gSD and (13:251) is found from (Oil . 

We will now briefly describe what we have referred to above as "less common" 
cases, depicted in Fig. |5l Let us first address the situation depicted in Fig. |5^, which 
occurs for, e.g., Az = 0.0044, 0.0045, 0.00474, 0.0049. The inner pair of instability 
peaks is found in these cases as described above. For example, for Az = 0.0049, it is 
obtained using {Np + Nm) = —2 and {Np — Nm) = 62. The corresponding calculated 
instability increment and the peak separation are PAp ^ 0.0051 and 2Q ^ 2 ■ 0.48, 
which are very close to the experimentally observed values. However, if for the same 
value of Az one takes {Np + Nm) = —3, then one finds that an instability can exist at 
Q ~ 0.78, corresponding to {Np — Nm) = 101. The calculated and observed values of 
the instability increment for this outer pair of peaks are 0.0030 and 0.0039, respectively. 

Let us note that we checked — by trial and error, as described in the previous para- 
graph — for the existence of a secondary pair of instability peaks for every value of Az 
reported in Figs. |3] and HI For each Az where such a secondary pair was observed in nu- 
merics, we also found it analytically, with the agreement between the calculated and ob- 
served values of Q being excellent and those of being good (the discrepancy between 
such values for Az = 0.0049, reported above, was the worst one we found). We even 
found, both numerically and analytically, small tertiary peaks for Az = 0.0044. On the 
other hand, we analytically found secondary peaks for Az = 0.0042, 0.00475, 0.00477 
where they were not observed in numerics at z = Zmax = 2000. However, at smaller 
z, such secondary peaks were indeed observed. The reason they were not observed 
for z = 2000 was that they were "drowned" by the pedestal of unstable modes which 
"rose" around the primary peaks due to modulational — i.e., non-numerical — insta- 
bility about the primary peaks. 

We will now comment on the calculation of the central peak for the case of Az = 
0.0050, shown in Fig. [5b. Here the value of (14. 7p with {Np + Nm) = —1 is approxi- 
mately 11.1. This is not negative, as we wanted this expression to be in the previous 
calculations (see the line after (14. 7p ). However, for {Np — Nm) = 0, the right-hand 
side of I K29\i is 11.1 G (4, 12), i.e. it is within the instability interval f H3|) with Q = 
(see (I3.28P ). According to the discussion found around Eq. (14. 3p . the instability peak 
a.t Q = should contain several grid nodes, which is confirmed by Fig. |5b. 

Finally, let us revisit the estimate for the order of magnitude of Q, which determines 
the frequency separation of the primary instability peaks. In (13. lip we stated, as an 
empirical observation, that Q = 0(1). We will now use Eq. (14. 8 p to generalize it. 
Note that for the largest integer value of {Np + Nm) that still renders expression (14. 7p 
negative, the expression in the square brackets in (14. 8 p is no greater than 2tt ■ 1/e. 
The corresponding {Np — A^j\/) estimated, and hence {Np — Nm), is then of the order 
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0{^/Tjt). Substitution of this into (13:28|) yields 

n = 0(l/v^). (4.9) 

Thus, as one makes the time window wider in the order-of-magnitude sense while keep- 
ing the other parameters fixed, the separation between the primary unstable Fourier 
modes should, on average, descrease. An implication of this for observing the numerical 
instability is presented in Section 5. 

5 Conclusions 

We reported, and then analytically explained, a numerical instability in the split-step 
Fourier method (11. 2p applied to the nonlinear Schrodinger equation (II. ip with the 
background solution being the soliton (II. 8p . Properties of this instability, such as the 
dependence of its increment and unstable mode's location on the step size, numerical 
domain's length T, and details of the background solution, are quite different from 
those properties on the background of the monochromatic wave (II. 4p . previously ob- 
tained by the von Neumann analysis in [6]. Namely, the dependence of the instability 
on those parameters is seemingly very irregular, as we illustrated in Tables 1, 2 and 
Figs, m m This is to be contrasted with monotonic dependences of the numerical insta- 
bility observed for finite-difference schemes in both constant- and variable-coefficients 
equations, as described in textbooks on numerical methods. Moreover, we demon- 
strated (see the end of Section 2) that the principle of "frozen coefficients" is not valid 
for the split-step Fourier method on the background of a localized nonlinear wave. In 
particular, the instabilities on the background of a single soliton, on one hand, and 
several well-separated identical solitons, on the other, can be drastically different. 

The analysis presented in Section 3 (see also Appendix 2) revealed that unstable 
modes can be found only near the resonant frequencies ±a;^, ±.uj2t,, etc. (see (II. 6p 
and the sentence after that equation). Interestingly, far away from the soliton, these 
unstable modes are not exactly periodic in the spatial variable t. Rather, their spatial 
envelope is exponentially growing or decaying in t (see Eqs. (13.80 . (I3.10p . (I3.15p . and 
(I3.22P ). What makes these modes satisfy the periodic boundary conditions, which are 
implicitly imposed by the use of Fourier transform in (II. 2p . is the change that they 
undergo near the soliton. Thus, the finite size of the numerical domain T is critical 
in our instability analysis, which is in stark contrast with the standard von Neumann 
analysis. 

As we noted above, the dependence of the instability increment on the parameters 
of the soliton and the numerical scheme is irregular, and there is no means to predict 
it "quickly", i.e., bypassing the procedure illustrated in Section 4.2. Again, this is in 
contrast with the instability analysis on a monochromatic-wave background [6] (see 
[8] and Appendix 3). However, one can still obtain two general conclusions: (i) 
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what maximum value this increment can have and (ii) how the chances to observe a 
numerical instability depend on the size T of the numerical domain. 

The maximum instability increment is obtained from Eqs. fl3.22p and fl3.25p by 
setting 2J'[[/2](o) + i{/3/-f){Jp + Jm) = 0. This yields: 

Interestingly, the last expression above is the same as an analogous expression in the 
monochromatic-background case (see Eq. (1A3.5I) in Appendix 3), where U{t) = A. 

Now let us show that as the time window length T is substantially decreased, the 
chances to observe numerical instability in any given simulation using fll.2p with f l2.ip . 
also decreasj§. This conclusion follows from a combination of arguments that led 
to formulae (14. 3p and f l4.9p . Indeed, recall from a discussion after Eq. f l4.8p that an 
instabihty would arise only if near the non-integer number {Np — A^m) estimated, there is 
an integer number {Np — N]\j) that would make the right-hand side of f l3.29p fit within 
the interval (14. 5p . A sufficient condition that would guarantee that such an {Np — Nm) 
can be found is obtained similarly to (14. 3p : 

47rfiestimated < 2 | (7//?) J^ft/^] (2nestimated) | , (5.2) 

where l^estimatod = {N p — N m) estimated / T . As follows from the discussion around (14. 9p . 
^estimated = 0{l/\/eT). Thus, as T decreases, the chances that condition (15. 2 p may be 
satisfied, also decrease. Hence the smaller T is, the "less often" a numerical instability 
of scheme (II. 2p on the background of a soliton would be observed. On the other hand, 
the smaller T is, the stronger the numerical instability, if it is observed, is on average; 
this follows from (I5.ip . Both these conclusions agree with our observations in Section 
2: Compare Table 1, obtained for T = 327r, with Fig. [3l obtained for T = 1287r. 
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Appendix 1: Derivation of Eq. ( 13.71 ) 



We verified that one can derive (13. 7p from (13.61) by using a Taylor series expansion of 
exp[i/3(ci;^ — uDAz] and (w„+i —Vn) in powers of e. However, an alternative derivation 
presented below is much less tedious and, importantly, more intuitive. 



^Note the word 'substantially'. As illustrated by Table 2 in Section 2, changing T by only a fraction 
of its value affects the occurrence of the instability in a non-monotonic and seemingly irregular way. 
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Let us first note that the last equation in (II. 2p is equivalent to the nonlinear 
Schrodinger equation (ll.ip plus a term proportional to 

Az/3j ■ [du, \u{t,z)\^]u{t,z) + 0{Az^), (Al.l) 

where [. . . , . . .] denotes a commutator. This follows from the Baker-Campbell-Hausdorff 
formula; see, e.g., Sec. 2.4.1 in [T]. Next, Eq. (13. 3 p is a linearized version of the 
last equation in (11.21) . Therefore, it must be equivalent to the linearized nonlinear 
Schrodinger equation plus terms of order 0{(3Azdtt) = 0{f3ijj'^Az) = O(e^), provided 
that we assume that w ~ 9^ = 0(1), or, equivalently, that the central frequency in 
expanding exp[il3uj'^Az] in a Taylor series is 0. 

In writing (13. 4p and then (13.60 . we stated that the central frequency is u-,, (or —u.,,) 
rather than 0. Correspondingly, Eq. (13. 6p must be equivalent to a modified linearized 
nonlinear Schrodinger equation written for a small deviation v, plus terms of order 
0(/3(aj^ — ul)Az) = 0(e) (see the text after (13. 4p ). Here the modification consists 
in replacing the operator du, whose Fourier symbol is — u;^ = — (u;^ ~ 0^), with the 
operator du + ool, whose Fourier symbol is — (w^ Thus, (13. 6p in the time domain 

is 

= -iPivtt + oolv) + hiul^v* + 2|usoip5) + 0(e) . (Al.2) 

Substituting into this equation Usoi from (II. 8p . changing the variable v = w exp[iKz], 
and neglecting the 0(e) term, one obtains Eq. (13.70 . 

Appendix 2: Location of instability peaks 

Here we will present an explaination of why the frequencies of unstable modes must 
be near ±uJt^, ±uj2n, etc. 

Suppose we seek the instability near a pair of frequencies ±a;o; i.e., we assume that 

\uj - cuol = 0{1) or |w- (-Wo) I =0(1). (A2.1) 

Then, proceeding as explained in the text after Eq. (13.30 . we obtain an equation similar 
to ( 13. 4p . where the "(—1)" and on the right-hand side are replaced with exp[i0o] 
and Uq, respectively, where 

00 = (3colAz . (A2.2) 
Then (13. 5p and (13.60 get replaced with 

^n = e-*<^«"u„, (A2.3) 

= e^^(-'--o)^^ T [vn + tlAziul,v: ■ e-^-^»" + 2\u,oi\'dn)] , (A2.4) 

where the term making the key difference between (13.60 and (]A2.40 is underlined. Let 
us now note that if the phase rotation, — 20o, in that term would equal —2'kN, 
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where N is any integer, then that term would equal 1, and the subsequent analysis 
would proceed as in Section 3.2 without any changes. Therefore, we can say that the 
nontrivial phase rotation in (1A2.4I) is — (20o — ^ttNq), where A^^o is the nearest integer 
to (Po/tt- For simplicity, but without loss of generality, let us assume that A^^o = 1; the 
case of Ai'o 7^ 1 is completely analogous. Then, Eq. (IA2.4P becomes 



= e^^(-'--o)^^ J- [v^ + ^lAz{ul,v: e'^^'^K-^)"^^ + 2|n,oir^„) J , (A2.5) 

where in rewriting the exponential term we have used (IA2.2P and (11. 6p . 

To go from the discrete equation flA2.5P to a counterpart of the continuous equation 



(13. 7p . we make two observations. First, nAz = z in the second exponential term in 
( IA2.5p . Second, and perhaps counter-intuitively: Despite the presence of this possibly 
fast-oscillating exponential, Eq. ( ]A2.5I) still describes a small change for {vn+i — Vn)- 
This is due to the presence of the small terms i'-fAz and f3{u)'^ — uJq)Az (see f lA2.1|) ) 
on the right-hand side of that equation. Therefore, the continuous variable v{t, z) 
interpolating the discrete variable in (IA2.5I) satisfies a counterpart of ( 1A1.2I) : 

Vz = -i[5{vu + ujIv) + ii{ul,v* e-2^^("o-"')^ + 2|nsoipi)) + 0(e) . (A2.6) 
In analogy with the argument presented in Appendix 1, we make a change of variables 

v = we^p [i{K - (5{ujI - ujI))z\ , (A2.7) 
which transforms (1A2.6I) into the following counterpart of (13. 7^ : 

Wz = -il3{wu + colw) -i{K- (5{ujI - ujI))w + i-^U'^iw* + 2w) . (A2.8) 



Then, a substitution analogous to (13. 8p with being replaced by ujq into Eq. ( IA2.8P 



yields a system of equations that is similar to (13. 9p . with the only changes being the 
replacements: 

uj^ by a;o and K by {K - P{ujI ~ ujD) . (A2.9) 

We now need to consider two cases: (i) P{oJq — = 0(1) and (ii) |/3(wo ~ 
col) \ ^ 1. We will show that in the first case, the results of analysis of Eqs. (13. 9 p with 
replacements ( ]A2.9P reduce to those obtained in Section 3.2, and in the second case, 



no instability can arise. 

In case (i), ujq — uJt, = 0(e), i.e. this case differs from that considered in Section 3.2 
only by a slight shift of the central frequency. Intuitively, such a shift could not change 
the location of the unstable peaks which we found to be away from by an amount 
of approximately Q = 0(1). Formally, this can be justified by a tedious calculation 
that reveals that Eqs. (I3.28p . (I3.29P with replacements ( 1A2.9P yield the same fl as the 
original Eqs. flX^ . fICT]) . Then, Eqs. f lXTU]) . ( KWf with replacements (^MM yield 
the same t-dependence of the solution of (13. 8p as the original Eqs. (I3.10p . (I3.15p . Thus, 
in case (i), the parameters of the instability reduce to those found in Section 3.2. 
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In case (ii), one cannot proceed as in case (i) by merely using replacements ( 1A2.9I) 
in Eqs. (13. 9p . The reason is that \K — (3{uJq — ^ 1, whereas in case (i) one had 
{K - f3{ui - uD) = 0(1). Indeed, in case (ii), a substitution (I3.10p with Q = 0(1) 
(which is our starting assumption — see (1A2.1I) ) would not yield psiow and msiow that 
would be slow functions of t; see the first term on the right-hand side of (13.141) . The 
only way the large term (^K — /3{uq — uj^)) could be eliminated from the couterpart of 
( 13. 9p is by using different ^-dependences in the exponentials in ( 13. 10^ : 

P = Psiow(^, z) exp [-int + 2i{/3/e)nz + /3Az - i{K - (5{ujI - ujI))z\ , (A2.10a) 

m = msiow(r, z) exp [int + 2i{l3/e)VLz + /3Az + i{K - /3(wo - ujI))z] . (A2.10b) 

(Note that the last terms in the exponents in ( 1A2.10I) essentially undo transformation 
(1A2.7I) .) In (1A2.10I) . Ps\ow and msiow are slow functions of t, but not of z. Substituting 
flA2.10l) into the counterpart of (13. 9p one would obtain, instead of the ^-independent 
system ()3.14p . a 2;-dependent system of the form: 

(Pslow)^ = fcl(Pslow)r + &2Pslow + h m^low e^* ('^"^^'^0"'^' ^) ^ (A2.11a) 

("^slow)2 = ci(msio„)^ + C2msiow + csPsiow e~^*('^"'^*^'^«"'^^^)^ (A2.11b) 

where all the coefficients bi through C3 are of order 0(1) and independent of z. The 
presence of rapidly oscillating exponential terms in ( ]A2.1ip makes the effect of the 
coupling terms negligible, and then system ( ]A2.1ip gets essentially decoupled into two 
independent equations for psiow and msiow, which does not exhibit any instability. Thus, 
in case (ii) numerical instability does not occur. 



Appendix 3: Instability on the background of a 
monochromatic wave 

Here we will use the method presented in Section 3.2 to find the location and growth 
rate of the numerically unstable Fourier modes of method (II. 2p on the background 
of a monochromatic wave (II. 4p with Qcw = 0. Let us note that these results can be 
obtained from formulae (65), (37), and (64) of [6J by expanding them in a power series 
of the step size Az (denoted there by r). In such a way, the growth rate of the most 
unstable mode was obtained in [8j. 

The starting point of our derivation is system (13. 9p . which holds both for the soliton 
and monochromatic- wave backgrounds. In the latter case, U{t) = A and K = 7^^, 
where without loss of generality we assume that A is real. Thus, now this system, 
unlike (13. 9p on the soliton background, has all constant coefficients, and hence we can 
look for its solution in the form 

{p, m} = {P, M} e'^''+^' , W = 0(1). (A3.1) 
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Note that unlike in (I3.10p . here the p- and m-components of the small deviation w 
have the same t-dependence. Also, we have used the notation W, not Q, in (lA3.ip . 



because, unlike Q, the variable W does not have the dimension of frequency. Rather, 
eW = W/u-n has the same dimension as Vt. 

Substitution of flA3.1|) into (13.91) with the aforementioned values of U and K yields: 



{2il3W - i-fA'^ + X)P- ijA^ M = 0, (A3.2a) 

- i-fA^ P + {2i/3W - i^A^ - A) M = 0, (A3.2b) 
where we have neglected terms O(e^). Then the instability growth rate is 



A = ^(7^2)2 - {2PW - 7^2)2 . (A3.3) 

The location of the unstable mode(s) follows from the definition of the mode's fre- 
quency, u; = — eW = Uj, — (W/uJt,), and the condition that the expression under the 
radical in (IA3.3P is positive: 

< /3iy < 7v4^ (A3.4) 
The maximum value of the growth rate occurs at the midpoint of this interval and is 

Amax = 7^'- (A3.5) 

Note that the periodicity condition for p and m does not play here a critical role 
in determining the instability increment, in stark contrast to the case of the soliton 
background considered in Section 3.2. Namely, as follows from (13.81) and ( lA3.ip . here 
this condition simply requires that the frequency u = — eW fall onto the frequency 
grid: — eW = 2ni/T, where i is an integer. Thus, if the width of the instability 
band is less than the frequency grid spacing: 

^A^ilPM < 27r/T, (A3.6) 

it is possible that 27r£/T may fall outside the instability band. In this case, the insta- 
bility will not occur even if Az exceeds the threshold (11.71) . This was originally pointed 
out by Weideman and Herbst [6] and studied in detail by Yang in [8] . 
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